
#include<iostream>
using namespace std;
#include<cmath>

/*
 * Newton 多项式, 用n阶多项式插值，通过(x_data,y_data)中的前 n+1 个数据点，返还 x 处的多项式插值 Pn(x)
 * Possible loopholes: x_data, y_data segmentation error
 * 			x_data[] have too dense values (<10^{-10}), can lead to large round-off error
 */
double NewtonPoly(int N, double *x_data, double *y_data, int n, double x){

	if( n== 0 ) return y_data[0];
	else{
		double a = y_data[n] - NewtonPoly(N, x_data, y_data, n-1, x_data[n]);
		for(int i=0;i<n;i++)
			a /= ( x_data[n] - x_data[i] );
		for(int i=0;i<n;i++)
			a *= ( x- x_data[i] );
		return NewtonPoly(N, x_data, y_data, n-1, x) + a;
	}
}

int main(){

	int n;
	double x,y,dy;

	double *x_data = new double [n];
	double *y_data = new double [n];

	FILE *fp=fopen("NewtonPoly_input.txt","r");
	if(fp==NULL){
		cout<<"error: failed to open NewtonPoly_input.txt.\n";
		exit(1);
	}

	n=0;
	while( fscanf(fp, "%*c%*[^\n]") != EOF ){
		n++;
	}
	cout<<"n="<<n<<endl;
	
	fp=fopen("NewtonPoly_input.txt","r");
	for(int i=0;i<n;i++){
		fscanf(fp, "%lf %lf", &x_data[i], &y_data[i]);
	}
	fclose(fp);

	double xmin = x_data[0], xmax = x_data[0];
	for(int i=0;i<n;i++){
		if( xmin > x_data[i] ) xmin = x_data[i];
		if( xmax < x_data[i] ) xmax = x_data[i];
	}
	xmin -= 2;
	xmax += 2;

	int num_grid = 100;
	double step = (xmax - xmin)/num_grid;

	fp=fopen("curve.txt","w");
	if(fp==NULL){
		cout<<"error: failed to open curve.txt.\n";
		exit(1);
	}

	for(int i=0;i<num_grid;i++){
		x = xmin + step*i;
		y = NewtonPoly(n, x_data, y_data, n-1, x);
		fprintf(fp,"%lf   %lf \n", x, y);
	}
	fclose(fp);

	return 0;
}

